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Abstract. While most chemical reactions in the interstellar medium take place in the gas phase, those occurring 
on the surfaces of dust grains play an essential role. Such surface reactions include the catalytic production 
of molecular hydrogen as well as more complex reaction networks producing ice mantles and various organic 
molecules. Chemical models based on rate equations including both gas phase and grain surface reactions have 
been used in order to simulate the formation of chemical complexity in interstellar clouds. For reactions in the 
gas phase and on large grains, rate equations, which are highly efficient to simulate, are an ideal tool. However, 
for small grains under low flux, the typical number of atoms or molecules of certain reactive species on a grain 
may go down to order one or less. In this case the discrete nature of the populations of reactive species as well 
as the fluctuations become dominant, thus the mean-field approximation on which the rate equations are based 
does not apply. Recently, a master equation approach, that provides a good description of chemical reactions 
on interstellar dust grains, was proposed. Here we present a related approach based on moment equations that 
can be obtained from the master equation. These equations describe the time evolution of the moments of the 
distribution of the population of the various chemical species on the grain. An advantage of this approach is the 
fact that the production rates of molecular species are expressed directly in terms of these moments. Here we use 
the moment equations to calculate the rate of molecular hydrogen formation on small grains. It is shown that the 
moment equation approach is efficient in this case in which only a single reactive specie is involved. The set of 
equations for the case of two species is presented and the difficulties in implementing this approach for complex 
reaction networks involving multiple species are discussed. 
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i-h ! 1. Introduction 



The chemistry of interstellar clouds consists of reac- 
tions taking place in the gas phase as well as on th e 
surfaces of dust grains flHartquist fe Williams 1995| ). 
Reactions that take place on dust grain surfaces in- 
clude the formation of molecular hydrogen as well as 
reaction networks producing ice mantles and various 
organic molecules. The formation of H2 on dust grain 
surfaces is a process of fundamental importance due 
to the fact that H2 cannot form in the gas phase effi- 
ciently enough to account for its observed abundance 



process of H2 for mation has attr act ed the attention of 
both theorists (Williams 1968; Bmoluchowski 1981 



Smoluchowski 1983| 


Aronowitz & Chang 1985 


Duley & Williams 1986^ 


Pirroncllo & Avera 1988 


Sandford & Allamandola 1993; Takahashi et al. 1999 


and experimentalists (Brackmann 196! 


; Schuttel976 


Pirroncllo ct al. 1997a 


Pirroncllo et al. 1997b 


Pirroncllo et al. 1999; 


Manico et al. 2001 


• 



(Gould fc Salpeter 1963; Hollenbach fc Salpeter 1970 



Hollcnbach fc Salpeter 1971; Hollenbach et al. 1971) 



The computational modeling of chemical reac- 
tion networks on dust grains in the interstellar 
medium is typically done using rate equation mod- 
els (|Pickles fc Williams 1977|; Id'Hendecourt 1985 ; 



Brown 1990 



Molecular hydrogen plays a crucial role in gravitational 
collapse and star formation by absorbing energetic 
photons and radiating the energy in the infrared, to 
which the cloud is more transparent, thus cooling it and 
enabling further collapse. Furthermore, H2 molecules 
are a necessary component for the initiation of chemical 
reaction networks that give rise to the chemical com- 
plexity observed in interstellar clouds. Therefore, the 



Hasegawa et al. 1992 



Brown fc Charnley 1990 



Hasegawa & Herbst 1993a; 



Hasegawa & Herbst 19931 


>; Casclli ct al. 1993 


Willacy & Williams 1993; 


Bhalabiea & Greenberg 1994) 



These models consist of coupled ordinary differential 
equations that provide the time derivatives of the densi- 
ties of the species involved. Integration of these equations 
provides the time evolution of the densities. As long as 
the grains are not too small and the flux is not too low, 
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the mean-field approximation applies. In this case, rate 
equations are an ideal tool for the simulation of surface 
reactions, due to their high computational efficiency. 
However, in the limit of very small grains under very 
low flux, rate equations are not always valid. This is 
because they take into account only average densities 
and ignore the fluctuations as well as the discrete nature 
of the populations of the atomic and molecular species 



(Tielens 1995 


Gharnley et al. 1997 




Caselli et al. 1998; 


|Shalabiea et al. 1998 




Stantcheva et al. 2001 


). These 



features become significant in the limit of small grains 
and low incoming flux, typically encountered in diffuse 
interstellar clouds where hydrogen recombination as well 
as various other grain surface reactions take place. For 
example, as the number of H atoms on a grain fluctuates 
in the range of 0, 1 or 2, the H 2 formation rate cannot 
be obtained from the average number alone. This can 
be easily understood, since the recombination process 
requires at least two H atoms simultaneously on the 
surface. 

Recently, a master equation approach was proposed, 
that is suitable for the simulation of chemical reactions on 



microscopic grains ( 


Biham et al. 2001 


Green et al. 2001 


Biham & Lipshtat 2002 


). It takes into account both the 



discrete nature of the H atoms as well as the fluc- 
tuations. In the case of hydrogen recombination, its 
dynamical variables are the probabilities -Ph(Ah) that 
there are N-r atoms on the grain at time t. The time 
derivatives -Ph(Ah), iVn = 0,1,2,... are expressed in 
terms of the adsorption, reaction and desorption terms. 
The master equation provides the time evolution of 
-Ph(Ah), Njj — 0,1,2,..., from which the recombina- 
tion rate can be calculated. The master equation ap- 
proach has been applied to the study of reaction net 



works involving multiple s pecies ( [Stantcheva et al. 2002 



Stantcheva fc Herbst 2002] ) . Since the state of the system 
is given by the set of number densities of all species, the 
number of equations quickly increases with the number 
of reactive species. Therefore, suitable cutoffs should be 
imposed in order to keep the number of equations at a 
tractable level. A Monte Carlo approach based on the 



master equation was also proposed (Charnley 2001), and 



applied for reaction networks involving multiple species 



(Stantcheva fc Herbst 2002) 



In this paper we present a set of moment equations 
that can be derived from the master equation, and use 
it to calculate of the rate of molecular hydrogen forma- 
tion on small grains. In the moment equations, the time 
derivatives of the moments (N^), k = 1,2,..., of the dis- 
tribution Ph(Nh) are expressed in terms of (Ah), (-^h), 
. . ., (A^ +1 ). To simulate the moment equations we impose 
a cutoff at a suitable value of fc, giving rise to k coupled 
linear differential equations. With such cutoffs, the mo- 
ment equations can be used efficiently to calculate the for- 
mation rate of molecular hydrogen for both steady state 
and time dependent conditions. Furthermore, they pro- 
vide a useful asymptotic expression for the rate of molec- 
ular hydrogen formation in the limit of very small grains. 



This expression can be used in order to integrate the H 2 
production rate over the grain-size distribution in an in- 



tcrstcllar cloud (Mathis 1977; Drainc 1985; Draine 1985 



Mathis 1990; Mathis 1996; Weingartner 2001; Cox 2002), 



and obtain the total production rate per unit volume of 
the cloud. 

In the case of complex reaction networks involv- 
ing multiple species, the number of moment equations 
quickly increases as the number of reactive species grows. 
Furthermore, the higher moments of the distribution are 
large and cannot be neglected. Thus, setting up the cut- 
offs in order to limit the number of equations is compli- 
cated and tends to introduce significant errors. As a result, 
for the case of multiple species we have not been able to 
develop an effective computational framework based on 
the moment equations. Here we present the set of mo- 
ment equations for two species, such as hydrogen and oxy- 
gen, where the moments take the form (N-^Nq), k,l = 
0,1,2,..., and discuss the approximations that may be 
used in order to set suitable cutoffs in these equations. In 
spite of these difficulties, the moment equation approach 
has some advantages that justify further attempts to de- 
velop it for the case of multiple species. One advantage is 
that the moment equations are a direct generalization of 
the rate equations and resemble their structure. One may 
thus expect that they will be suitable for incorporation 
into models of interstellar chemistry based on rate equa- 
tions. Another advantage is that in this approach, the mo- 
ments, which quantify the production rate of molecules, 
are directly calculated unlike the case of the master equa- 
tion where the calculation of moments requires further 
processing. 

The paper is organized as follows. In Sect. 2 we con- 
sider reactions involving a single atomic specie, focusing 
on the case of hydrogen recombination. In Sect. 3 we con- 
sider reactions involving two atomic species, such as hy- 
drogen and oxygen. In each of these Sections we briefly 
describe the rate equations and master equation for the 
system under study and then introduce and analyze the 
moment equations. The results are discussed and summa- 
rized in Sect. 4. 

2. Hydrogen recombination 

Consider a diffuse interstellar cloud dominated by a den- 
sity ph (cm -3 ) of H atoms, and includes some density of 
dust grains. The typical velocity wh (cm s _1 ) of H atoms 
in the gas phase is given by 



8 fcs7g as 
7r ran 



(1) 



where m# = 1.67 • 10 24 (gram) is the mass of an H atom 
and T gas is the gas temperature (Landau & Lifshitz 1980). 



To evaluate the flux of atoms onto grain surfaces we will 
assume, for simplicity, that the grains are spherical with 
a radius r (cm). The cross section of a grain is a = nr 2 
and its surface area is 4o\ The number of adsorption sites 
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on a grain is denoted by S. The flux Fh (atoms s" 1 ) of 
H atoms onto the surface of a single grain is given by 
Fr — ph^h 17 - This flux can also be expressed according 
to F H = / H • S, where / H (ML s -1 ) is given by / H = 
/°h u h/(4s), and s (sites cm -2 ) is the density of adsorption 
sites on the surface. The H atoms stick to the surface and 
hop as random walkers between adjacent sites until they 
either desorb as atoms or recombine into molecules. The 
desorption rate of an H atom on the surface is 



W u = vexp(-E 1 /k B T), 



(2) 



where v is the attempt rate (typically assumed to be 10 12 
s -1 ), Ei is the activation energy barrier for desorption of 
an H atom and T is the surface temperature. The hopping 
rate of an H atoms is 



a H 



exp(-F /fc B T), 



(3) 



where Eq is the activation energy barrier for H diffusion. 
Here we assume that diffusion occurs only by thermal 
hopping, in agreement with recent experimental results 
(Katz et al. 1999). Throughout this paper we use the pa- 
rameters obtained experimentally for amorphous carbon, 
namely the activation energies are Eq = 44.0 meV and 
E\ = 56.7 meV (Katz et al. 1999), and the density of ad- 
sorption sites on the surface is s ~ 5 x 10 13 (sites cm~ 2 ) 



(Biham et al. 2001). For the density of hydrogen atoms 
in the gas phase we take pn = 10 (atoms cm -3 ). The gas 
temperature is taken as T gas = 90 K, thus vh — 1-37 x 10 5 
(cm s -1 ). These are typical values for diffuse interstellar 
clouds. 

The number of H atoms on the grain is denoted by Ah 
and its expectation value under the given conditions is de- 
noted by (-/Vh). The rate An = clh/S is approximately the 
inverse of the time t s required for an atom to visit nearly 
all the adsorption sites on the grain surface. This is due 
to the fact that in two dimensions the number of distinct 
sites visited by a random walker is linearly proportional 
to the number of steps, up to a logarithmic correction 
QMontroll fc Weiss 19651) . 



2.1. Rate equation 

The rate equation describing H2 formation on dust grain 
surfaces takes the form 



d(N H ) 
dt 



Fh-W h (N h )-2A h (N h ) 



(4) 



The first term on the right hand side describes the flux 
of H atoms, the second term describes the desorption of 
H atoms and the third term describes the diffusion and 
recombination. Here we assume, for simplicity, that all H2 
molecules dosorb from the surface upon formation. The 
production rate i?^ ram (molecules s -1 ) of H2 molecules 
from a single grain is given by 



As long as (-/Vh) S> 1 Eq. (||) provides a very good evalu- 
ation of the H2 formation rate. However, for very small 
grains and low flux, for which (-/Vh) — 1 or less, this 
mean-field approach fails, and tends to overestimate the 
H 2 production. This is because the rate equation does not 
take into account the discrete nature of the population of 
H atoms, and the fact that it takes at least two atoms 
on the surface simultaneously to enable recombination. In 
this limit, the master equation approach is needed in order 
to evaluate the H2 formation rate. 

2.2. Master equation 

We will now describe the master equation for H2 formation 
on small interstellar grains, exposed to a flux Fh of H 
atoms. At any given time the number of H atoms adsorbed 
on the grain may be ./Vh = 0, 1, 2, . . ., and the probability 
that there are ./Vh hydrogen atoms on a grain is given by 
Fh(Ah), where 



J2 MNh) = 1. 



(6) 



«H=0 



The time derivatives of these probabilities, Ph(Ah), are 
given by 



Pn(l) 



Fh(-ZVh) 



F|f n = A h (Nh) 2 . 



(5) 



-FhFh(0) + W h Fh(1) + 2 ■ 1 • A h Fh(2) 
Fh[F h (0)-F h (1)] 

Wn [2F H (2) - F H (1)] + 3 • 2 ■ A H F H (3) 

Fh [Fh(-/Vh — 1) — Fh(-/Vh)] 

W B [(./Vh + 1)P H (A H + 1) - -/Vh-Ph(-ZVh)] 

Ah [(-/Vh + 2)(A H + 1)F H (/V H + 2) 

A H (A H - l)fk(JV H )]. (7) 



Each of these equations includes three terms. The first 
term describes the effect of the incoming flux Fh on the 
probabilities. The probability Ph(-ZVh) increases when an 
H atom is adsorbed on a grain that already has -/Vh — 1 ad- 
sorbed atoms [at a rate of FhPh(Ah — 1)L and decreases 
when a new atom is adsorbed on a grain with ./Vh atoms 
on it [at a rate of FhPh(-ZVh)]- The second term includes 
the effect of desorption. An H atom that is desorbed from 
a grain with Ah atoms, decreases the probability Ph(Ah) 
[at a rate of NnWnPn(Nn)], and increases the probabil- 
ity Ph(-ZVh — 1) at the same rate. The third term describes 
the effect of recombination on the number of adsrobed H 
atoms. The production of one molecule reduces this num- 
ber from /Vh to Ah — 2. For one pair of H atoms the recom- 
bination rate is proportional to the sweeping rate An mul- 
tiplied by 2 since both atoms are mobile simultaneously. 
This rate is multiplied by the number of possible pairs of 
atoms, namely Ah (Ah — 1)/2. Note that the equations for 
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Ph(0) and Ph(1) do not include all the terms, because at 
least one H atom is required for desorption to occur and 
at least two for recombination. The rate of formation of 
H 2 molecules, i?§ 2 am (molecules s -1 ), on the surface of a 
single grain is thus given by 

- (JV H >] 



i?| r 2 am = A 



(8) 



where 



(9) 



N H =0 



is the fcth moment of the distribution. 

The time dependence of BS™ can be obtained by nu- 
merically integrating Eqs. (J7]) using a standard Runge- 
Kutta stepper. For the case of steady state, namely 
Pn(Nn) = for all A?h, an analytical solution for Ph(Ah) 
is available, given in terms of Ah/Wh and Wr/Fr 
( preen et al. 2001| ; |Biham fc Lipshtat 2002| ). 



2.3. Moment equations 

By expanding the time derivatives (iVjj), k = 1, 2, . . ., us- 
ing Eq. (||) and inserting the expression for Ph(A h ) from 
the master equation (|t]), one obtains the moment equa- 
tions- 



tions: 

d(N K ) 
dt 

dt 
dt 



d(N& 
dt 



= F H + (-W H + 2A H )(N H ) - 2A K {N£) 

= F H +(2F H +W il -4A H )(N H ) 
+ (8A H - 2Wn){N n ) - 4A B (N&) 

= F H +(3F H -W U + 8A H )(N H ) 

+ (3P h + 3W h -2(L4h)<A£) 

+ (18A H - 3Wk)(iV|> - 6A H «) 

= F H ((l + 7V H ) fc -0 

+ ^ H (iVH[(^H-l) fe -<]) 

+ An (-/Vh (iVu - 1)[(JV H - 2) fc - iVg]>. 



(10) 



This is a set of coupled differential equations, that are 
linear in the moments (iVj|). The equation for the fcth 
moment depends on all the moments up the the (k + l)th 
order. A numerical integration of these equations would 
provide the values of these moments at any given time. 
The H2 production rate, that depends only on the first and 
second moments, could then be obtained. However, the 
difficulty is that we need to truncate this set of equations, 
say after the equation of the fcth moment. The (fc + l)th 
moment appears in this equation and we end up with more 
unknowns than equations. The (fc + l)th moment should 



then be approximated as a function of the first fc mo- 
ments. The approximation should be examined carefully 
since {N^ +1 } is not a small quantity (unlike the case when 
the master equation is truncated). In fact, due to the dis- 
creteness of TVh, the moments increase monotonically with 
their order, namely (iVg) < (A£ +1 ), fc = 1,2, . . . (and the 
< symbol is replaced by the < symbol if Ph(Ah) > for 
at least one value of Ah > 1). 

2.3.1. Setting the cutoffs 

Consider the case in which only up to fc hydrogen atoms 
are allowed to reside simultaneously on the surface of a 
grain. The moments of the resulting truncated distribution 
P H (A H ) will be 



(JVh) 

(JVh> 
<7Vh +1 > 



2P H (2) 
4P h (2) 



fcP H (fc) 
fc 2 P H (fc) 



2 fe P H (2) + ... + fc fe P H (fc) 

2 fe+1 P H (2) + ... + fc' £ + 1 PH(fc). (11) 



The first fc equations in Eqs. ([11]) can be used in order to 
express P H (1) ,. . ., P H (fc) in terms of (N£), . . ., The 
results can then be plugged into the (fc + l)th equation in 
Eqs. (pi]|), thus expressing (N^ +1 ) in terms of the first fc 
moments. For example, in the case of fc = 2 we obtain the 
third moment as a function of the first two: 



<A£) = 3<A£)-2<A£>. 



(12) 



Similarly, when Ph(Ah) is truncated after fc = 3 or fc = 4 
we obtain 

(N*) = 6(N&)-U(N& + 6{N*) (13) 
and 

(N H ) = -24(iVA) + 50(N H ) - 35<A|) + 10(N*), (14) 

respectively. In general, when up to fc hydrogen atoms 
are allowed to reside on a grain the (fc + l)th moment, 
(A^ +1 ), can be expressed as a linear combination of the 
first fc moments according to 



(N* +1 ) = '£c k+1 (n){NZ). 



(15) 



n=l 



The coefficients Ck+i(n) are obtained from the solution of 
a set of linear algebraic equations, that can be expressed 
in a matrix form as 



VC = v. 



(16) 



The matrix V is the Vandermonde matrix of size fc x fc, 
namely, V mn = m", where m, n = l,2,...,fc, and the 



Lipshtat & Biham: Moment equations for reactions on dust grains 



5 



vector v consists of 



m = 1,2,. 



,fc. The 



desired coefficients are the elements of the vector C = 
[Cfe + i(l), . . . , Cfc + i(fc)]. The coefficient Ck+i(n) turn out 
to be equal to the coefficient of x n in the polynomial 



(Bender et al. 2002) 



.(*) 



-n< 

3=0 



(17) 



Having an expression for (N^ +1 ) as a linear combination 
of the first k moments, we can now insert it into the last 
equation in Eq. (|l(]). We then obtain a set of k coupled 
linear differential equations for the first k moments. This 
set of equations can be simulated using a standard Runge- 
Kutta or any other stepper routine. A convenient choice of 
initial condition may be the case of an empty grain. In this 
case the moment equations are initiated with (N^) = 
for all n > 1 (while the master equation is initiated with 
P H (0) = 1 and P H (n) = for n > 1). 

2.3.2. Calculations and results 

To determine the desirable value of the cutoff k, for steady 
state calculations using the moment equations, it is useful 
to first calculate (iVk) using the rate equations [Eq. (Q)] 
for the same parameters used in the moment equations. 
We found that a suitable choice for the truncation is k = 
K-^h) + CI where \x] is the smallest integer which is 
larger than x. The parameter C is determined such that 
a sufficient number of equations are included to provide 
a good agreement with the master equation. In all the 
calculations presented here we used C = 1.2. 

We have used the moment equations to simulate the 
hydrogen recombination process on a grain under steady 
state conditions, and compared the results with those of 
the rate equations and the master equation. The surface 
parameters used are of amorphous carbon. The gas density 
and temperature are those specified above, and the surface 
temperature is T = 18 K. 

The number of H atoms on the grain is shown in Fig. 
0(a) as a function of the number of adsorption sites, S, 
using the rate equations (dashed line), master equation 
(O) an( i moment equations (solid line). The production 
rate of H2 molecules on the surface of a single grain is 
shown in Fig. |l](b) as a function of S. The results of the 
rate equations (dashed line) are simply linearly propor- 
tional to S. For large enough grains the results of the rate 
equations coincide with those of the master equation (O) 
and the moments equations (solid line). For small grains, 
the results of the moment equations are in perfect agree- 
ment with the master equation, while the rate equations 
greatly over-estimate the production rate. In the simula- 
tions of the moment equations, the cutoff was determined 
according the the procedure specified above. The number 
of equations used in each range of grain sizes is speci- 
fied. Note that using this procedure we obtain a perfect 
agreement between the moment equations and the mas- 
ter equation for the entire range of grain sizes. Results 
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Pi 

a 10 



Rate Equation 

Moment Equations 

Master Equation: 
O Complete set 
x Cutoff after 3 equations 
+ Cutoff after 5 equations 




O : 
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Fig. 1. (a) The expectation value of the number of H 
atoms on the grain as a function of the number of ad- 
sorption sites, S, on the grain, obtained from the moment 
equations, the master equation and the rate equations; 
(b) The production rate R^ ln of molecular hydrogen on 
a single grain of S surface adsorption sites. The results 
of the moment equations (solid line) are in perfect agree- 
ment with those of the master equation (O), while the 
rate equations (dashed line) over-estimate the production 
rate for small grains. The number of moment equations 
used for each range of grain sizes is specified. 



are also shown for the master equation with cutoffs af- 
ter the third equation (x) and after the fifth equation 
(+). It is shown that while the sets of two, three and four 
moment equations provide a perfect agreement with the 
complete master equation in the specified domains, the 
master equation with cutoffs (including at least as many 
equations) exhibits significant deviations. 

In the limit of very small grains, only two moment 
equations are needed in order to obtain perfect agreement 
with the master equation. These equation are derived from 
the first two moment equations in Eqs. ( |l0| ) and the cutoff 
condition given in Eq. (113). They take the form 
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Fig. 2b 

Fig. 2. (a) The expectation value for the number of H 
atoms on the grain vs. time; and (b) the production rate 
i?H 2 a °f m °l ecuiar hydrogen vs. time, for a large grain 
(S = 62,831 sites) and for a small grain (S = 628 sites). 
The initial condition is of an empty grain surface, while 
the surface parameters, the temperature and the flux are 
identical to those used in Fig. 1. For the large grain, it 
is found that the moment equations (solid line), master 
equation (O) and the rate equations (dashed line) are 
all in perfect agreement even during the transient time. 
For the small grain, the set of two moment equations is 
in perfect agreement with the complete master equation, 
while the rate equations deviate significantly. 



d{N K ) 
dt 

dt 



= Fu 



(2A H 
(2F H 



-W H ){N ll )-2A li (N 2 } 
W H + 4A H )<iVH) -(4A H 



vs. time is shown in Fig. 2(a), and the production rate 
i?H 2 am of molecular hydrogen is shown in Fig. 2(b). The 
initial condition is of an empty grain surface, while the sur- 
face parameters, the temperature and the flux are identical 
to those used in Fig. 1. For the large grain, it is found that 
the moment equations, master equation and the rate equa- 
tion are all in perfect agreement even during the transient 
time. For the small grain, the set of two moment equations 
is in perfect agreement with the complete master equation, 
while the rate equation deviates significantly. For steady 
state conditions, Eqs. (|l^) can be solved exactly, giving 
rise to 



(Nk) = 

<^H> = 



F H (A H + W H ) 



WhA h + W 2 
f Ah + Wh) 



2A H F H + WhAh 



(19) 



A simple exact expression for the production rate i?^™", 
of molecular hydrogen, given by Eq. (^J) , can now be ob- 
tained for steady state conditions in the limit of small 
grains and low flux. It takes the form 



R gram = 

n-2 



A 



S 2 



_ 1 + ^7Wh + 2 



(20) 



We observe that, as long as S is larger than both s v ; S i t = 
&h/Wh (the number of sites that an H atom typically 
visits before it desorbs) and s vacan t = Wh//h (the typical 
number of vacant sites per adsorbed H atom) , the produc- 
tion rate is linearly proportional to S, and the grain size 
does not play any special role ( Biham fc Lipshtat 2002 ). 



In the case of small grains for which S is smaller than 
both s v i s it and s vaC ant the production rate is reduced and 
becomes proportional to S 2 . The production efficiency 
V (S) = 2Rf ain 



t](S) 



2/h 



h 2 /Fa then takes the form 
S 



1 + 



s 



(21) 



We find that these equations are dynamically stable and 
can be used for time-dependent simulations. In Fig. 2 we 
present the results of time-dependent simulations using 
Eqs. ( |l8[ ) for a large grain (S — 62, 831 sites) and for a 
small grain (S = 628 sites). The expectation value (N-g) 



a H /W H ' Wh//h. 

In Fig. 3 we present the production efficiency, given by 
Eq. (|2"l|), as a function of S under the same conditions 
as those used in Fig. 1. For very small grains, namely 
for 5" < min{ Svisit, Svacant}, the efficiency r)(S) is linearly 
proportional to S, while for S > max{s v i s ;t, s vacant} it 
saturates and coincides with the rate equation result that 
is independent of S. 

These results can now be used in order to evaluate the 
total production of molecular hydrogen per unit volume 
of the cloud from grains of all sizes. Consider a cloud in 
which the grain size distribution is given by p S r(r) (cm -4 ) 

2WH)(li$$ xe ran S e rmin < r < r max, where r min and r max are the 
lower and upper cutoffs of the distribution [here p gr {r)dr 
(cm -3 ) is the number density of grains of sizes in the range 



(r,r+dr)]. The formation rate i?n 2 ( cm 
ular hydrogen in the cloud is given by 



Ra 2 — PhVh 



nr 2 p gI (r)r](S)dr. 



s ) of molec- 



(22) 
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entire range of grain sizes 



< r < rv, 



Unlike the 




Fig. 3. The efficiency tj(S) of molecular hydrogen for- 
mation on a grain of S assorption sites is presented. It 
is shown that in the limit of small grains the efficiency 
depends linearly on S, while for large grains it saturates 
towards the efficiency predicted by the rate equations and 
becomes independent of S. 




6x10 8x10 

r (cm) 



Fig. 4. The differential production rate of molecular hy- 
drogen, given by the integrand in Eq. ( ^2|) vs. the grain 
radius, r. The largest contribution comes from very small 
grains, of sizes around 20 nm. For larger grains the con- 
tribution decreases due to the fact that their total sur- 
face area is smaller. For smaller grains the contribution 
decreases due to the reduction in the recombination effi- 
ciency. 

where rj(S) is given by Eq. and S = 4vrr 2 s. 

To demonstrate these ideas we will consider a simple 
example in which the grain size distribution exhibits a 
power-law behavior of the form 



( \ K 



(23) 



where a = 3. The special property of this distribution 
is that the total grain mass is equally distributed in the 



mass, the total surface area of all grains in the size range 
between r and r + dr, scales like r _1 , namely most of the 
grain surface area is on the small grains. The prefactor 
K is determined by the condition that the total number 
density of grains is smaller than pn by a factor of 10 -12 . 
Under this condition K = 2 ■ 10~ n /(r~ 2 n _ r m 2 ax)- I n this 
case the integral can be easily evaluated and the result is 



i?H, = 



^Kp-av-af-RS ( Br, 



W H B 



In 



f 



Br 2 - 

mm 



f 



(24) 



where B = 4ns[Wii/aH + 2/h/Wh]- To evaluate the con- 
tribution of grains of different sizes to the total production 
rate of molecular hydrogen, we show in Fig. 4 the differ- 
ential production rate, given by the integrand in Eq. ( P2|) , 
vs. the grain radius, r, where r m ; n = 2 x 10~ 7 (cm) and 
r max = 1.25 x 10~ 5 (cm). Starting from the limit of large 
grains, the differential production increases with decreas- 
ing grain size due to the larger total surface area of smaller 
grains. However, below some grain size the recombination 
efficiency decreases and the differential production rate 
starts to decrease. We observe that under these condi- 
tions the largest contribution comes from extremely small 
grains of sizes around 20 nm. Observation indicate that 
the exponent in Eq. ( p3| ) is a — 3.5 (see e.g. Weingartner 
& Draine 2001). Therefore, the contribution of very small 
grains to molecular hydrogen formation may be even more 
dominant than in the analysis shown here. 

Observations indicate that the temperatures of dust 
grains in diffuse clouds are typically in the range be- 
tween 10 - 20 K. The dependence of the molecular hy- 
drogen formation process on the surface temperature of 
the grains is presented in Fig. 5. The average number 
of H atoms on the surface of a very small grain of ra- 
dius r = 10 -6 (cm) is shown in Fig. 5(a). The production 
rate of molecular hydrogen on such grain is shown in Fig. 
5(b). The results of the moment equations (solid line) and 
the master equation (O) coincide perfectly, while the rate 
equations (dashed line) significantly overestimate the pro- 
duction rate for temperatures higher than 16 K. As the 
temperature increases above T = 16 K, desorption be- 
comes faster and the production rate decreases. Below 16 
K the production rate saturates and becomes independent 
of the temperature within the temperature window of 12 
- 16 K. Below T = 12 K the production rate decreases 
sharply due to the lower mobility of H atoms on the sur- 
face flKatz et al. 1999|). 



3. Reactions involving two Species 

Observations in the past two decades provide evidence 
that chemical reaction networks on dust grain surfaces 
go much beyond the production of molecular hydrogen. 
In particular, ice mantles that consist of H2O and CO, 
as well as more complex organic molecules, form on grain 
surfaces. 

Here we consider a surface reaction network involving 



hydrogen and oxygen (Casclli et al. 1998). In this system 
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3.1. Rate equations 




The rate equations that describe the chemical reactions 
on dust grains in our simplified hydrogen-oxygen system 
are 



Fig.5a 



Moment Equation 

O Master Equation 
Rate Equation 




Fig.5b 



T(K) 



Fig. 5. The average number of H atoms on the grain 
surface (a), and the production rate (b) for a grain of 
radius r = 1CP 6 (cm) as a function of the grain tempera- 
ture. The results of the moment equations (solid line) and 
the master equation (O) coincide perfectly, while the rate 
equations (dashed line) significantly overestimate the pro- 
duction rate for temperatures higher than 16 K. As the 
temperature increases above T = 16 K, desorption be- 
comes faster and the production rate decreases. Below 16 
K the production rate saturates and becomes independent 
of the temperature within the temperature window of 12 
- 16 K. 



surface reactions produce H2, O2, OH and H2O molecules. 
While the hydrogen molecules typically desorb upon for- 
mation or shortly later, the heavier species tend to remain 
on the surface. Some of them, such as OH may then par- 
ticipate in further reactions. However, in the model pre- 
sented below we assume, for simplicity, that the reaction 
products do not participate in further chemical reactions 
(namely, H2O is not produced). In practice, we consider 
these molecules as if they all desorb from the surface upon 
formation. 



d(N H ) 
dt 

d(Np) 
dt 



Fr-W h (Nh)-2Ah(Nh) 2 
(A K +A ){N }i )(No) 
F -Wo{No)~2Ao(No) 2 
(A H + A )(N li )(No} 



(25) 



where Ao is the number of oxygen atoms on the grain 
and Ao = cio/S is their sweeping rate. Here we denote 
the density of the O atoms in the gas phase by po- Then- 
velocity, vo, is determined by Eq. (|l|) (replacing mjj by 
the mass mo of an oxygen atom) with the same temper- 
ature as that of the hydrogen gas. The flux of O atoms 
adsorbed on the grain surface is given by F<j = povoa 
(atoms s _1 ). The desorption rate Wo is given by an ex- 
pression analogous to Eq. (Q) with the same v and with 
Ei replaced by the activation energy for desorption of O 
atoms. The hopping rate ao of the O atoms is given by 
an expression analogous to Eq. (||), with Eq replaced by 
the activation energy for hopping of O atoms. 

As long as the grains are not too small and the fluxes 
are not too low, such that Ah 3> 1 and Ao 3> 1 these rate 
equations evaluate correctly the formation rates of all the 
reaction products. However, when the typical numbers of 
H or O atoms on the grain go down to order one, these 
equations become unsuitable and exhibit significant devi- 
ations from the actual production rates. 

3.2. Master equation 

The master equation that describes the hydrogen-oxygen 
system consists of a two dimensional matrix of equations 
for the time derivatives of Ph&o(Ah,Ao). It takes the 
form 



-Ph&o ( A H ,A ) =F h [ j Ph&o(Ah-1,Ao) 

- -Ph&o(A h ,A )] 

- Ph&o(A h ,A )] 

+ Wh[(A h + l)P H &o(A H + 1, A ) 

- A h Ph&o(A h ,Ao)] 

+ W [(N + l)P H &o(A H , Ao + 1) 

- A Ph&o(Ah,A )] (26) 
+ A H [(A H + 2)(A H + 1)Ph&o(A h + 2, A Q ) 

- Ah(A h -1)Ph&o(Ah,Ao)] 

+ A [(N + 2)(N + l)P H &o(A H , A + 2) 

- A (A -1)Ph&o(Ah,Ao)] 
+ (Ah + Aq) 
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[ (N H + 1)(N + l)P H &o(iVH + 1, N + 1) 
- N h NoPk&o(Nk,N )}, 

where Nn,No = 0,1,2,.... In the equations in which 
TVh = 0, 1 or iVo = 0, 1, some of the terms vanish. These 
terms can be easily identified since their initial or final 
states include a negative number of H or O atoms on the 
grain. 



3.3. Moment equations 

The moment equations for the hydrogen-oxygen system 
take the form 



d(N H ) 
dt 

d(N ) 
dt 

d(N£) 
dt 



dt 



= F H + (2A H - Wu)(N u ) - 2A H {N&) 

- (A H + A ){N H N ) 

= F Q + (2A Q - W ){N ) - 2A {N ) 

- (A li + A )(N li N ) 

= F H + (2F H + Wh - ^hX^Vh) 
+ (-2VK H + 8A H )(iV2} 

- 4A H (A^} + (A H + A )(N H N ) 

- 2(An + Ao)(NlN ) 

= F Q + (2F + W Q - lAo){N Q ) 
+ (-2W + SA ){N^) 

AAo(Nl) + (A H + A )(N H N ) 

- 2(A H + A )(N H N^ 

= F K (N ) + FoiNn) 

+ (-W H - Wo + 3A H + 3A )(N }i No) 
~ {3Ak + Ao)(N*No) 

- {A S +3A ){N H N&. (27) 



and the production rates of the molecular species are given 

by 



d(NuNc 
dt 



^gram 

ii2 
^grain 

U2 

ngrain 
n OH 



A H «2Vg) - (JV H » 
A {(Nl) - (N )) 
(A H + A )(N H N ). 



(28) 



In the moment equations the time derivative of each mo- 
ment is expressed as a linear combination of other mo- 
ments. The set can be extended to include higher mo- 
ments. However, the number of unknowns (namely mo- 
ments to be computed) grows faster than the number of 
equations. In Eqs. (j27j), consisting of five equations, there 
are nine unknowns: (N H ), (N ), (JVg), (iVg), (N&), (N Q ), 
{N H N ), (N&No) and (N H N^). 

Therefore, in order to solve the moment equations 
one has to use approximations for four of the unknowns 



listed above. In principle, these approximations can be ob- 
tained using the following properties. The number density 
of hydrogen in the gas phase is by orders of magnitude 
higher than that of any other specie, including oxygen. 
Therefore, Fh 3> Fq. The activation energies for diffusion 
and desorption of H atoms on the grain surface are much 
lower than those of heavier atomic species such as oxy- 
gen. Therefore, the diffusion of H atoms on the surface is 
much faster and dominates the chemical activity on the 
surface. The density of H atoms on the surface is thus pri- 
marily controlled by the processes of H2 formation and H 
desorption. The further reduction of H density due to the 
formation of OH can be considered as a small correction. 
We have examined such approximations and found that 
they do not work very well in practice and that the results 
deviate significantly from those of the master equation. 

4. Summary 

We have introduced a set of moment equations for the 
analysis of the formation of molecular hydrogen and other 
chemical reactions on dust grain surfaces in the interstel- 
lar medium. These equations are derived from the master 
equation that describes these processes. Like the master 
equation, the moment equations are exact even in the limit 
of small grains and low flux where the rate equations fail. 
They take into account the discrete nature of the popula- 
tions of atomic and molecular species that participate in 
the reactions. While the master equation is expressed in 
terms of distributions such as Ph(Nh), Nn = 0, 1, 2 . . ., 
the moment equations are expressed in terms of the mo- 
ments of these distributions. 

We have shown that the moment equation approach is 
useful for simulations of molecular hydrogen formation. 
We expect that this approach will make it possible to 
efficiently calculate the total production rate of molecu- 
lar hydrogen in diffuse clouds. The results will then be 
used in order to determine whether the molecular hydro- 
gen formation process studied here is efficient enough to 
account for the observed abundance of H2 . The input re- 
quired for such calculations includes the surface param- 
eters (obtained from laboratory experiments), the grain 
size distribution, the density and temperature of the H 
gas, the grain temperature as well as the dissociation rate 
of H2 molecules in the gas phase. Since the formation rate 
of molecular hydrogen is highly sensitive to the grain tem- 
perature and size distribution, the accuracy of the results 
may be limited by the quality of the observational data. 

For more complex reaction networks involving multi- 
ple species, setting the cutoffs in the moment equations 
requires certain approximations that introduce significant 
errors. Therefore, for system that consist of more than 
a single reactive specie, we currently do not have a use- 
ful computational framework based on the moment equa- 
tions. Potential advantages of the moment equations, jus- 
tifying further attempts to further develop this approach, 
include the fact that they are similar in structure to 
the rate equations, possibly making it easier to incorpo- 
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rate them into models of interstellar chemistry. Moreover, 
they directly provide the first and second moments, which 
quantify the production rates of the molecular species in- 
volved. This is unlike the master equation in which the 
moments are calculated from the distribution in a post- 
processing stage. 
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